Import modules



In [17]:

    
import sympy
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
from sympy import init_printing
from sympy.utilities.lambdify import lambdify
init_printing()
from JSAnimation.IPython_display import display_animation
from matplotlib import animation

Important constants



In [18]:

    
V_max = 136 # km/hr. The maximum velocity of cars on the road.
L = 11 # km. The length of the stretch of road.
rho_max = 250 # cars/km. Maximum density of cars before traffic stops.
nx = 51 # number of grid points in the x direction
dt = 0.001 # hours.
sigma = 0.5 # CFL number we use.
BG_VAL = 20 # The background value used for initial data (see below)

Useful functions



In [19]:

    
hours = lambda minutes: minutes/60.
minutes = lambda hours: hours*60.
mps = lambda kph: (1000./3600.)*kph # cast km/h to m/s
V = lambda rho: V_max * (1 - (rho/rho_max))
flux = lambda rho: V(rho)*rho

Grid



In [20]:

    
x = numpy.linspace(0,L,nx)
dx = x[1]-x[0]
print "dx = ",dx









    



dx =  0.22

Integration Methods



In [21]:

    
def rhs(rho):
    """
    given the equation
    d rho/dt = f(rho)
    returns f(rho)
    
    Inputs:
    --- rho, a vector of floats of lenth nx. 
    
    Returns:
    drho/dt on the grid excluding the zeroth element, which is not set
    """
    F = flux(rho)
    return -(F[1:]-F[:-1])/dx



In [22]:

    
def time_step(rho,rho0):
    """
    Given rho at time t, return rho at time t+dt
    
    Inputs:
    --- rho, rho(x) at time t
    --- rho0, the boundary condition, rho(x=0) for all t
    """
    out = numpy.empty_like(rho)
    out[1:] = rho[1:] + dt*rhs(rho)
    out[0] = rho0
    return out



In [23]:

    
def initial_data(background):
    """
    Given an array x, returns an array rho(x) at time 0. 
    Initial configuration is a tophat function
    
    rho(x[i]) = 50 if 10 <= i <= 20
    rho(x[i]) = background otherwise
    """
    rho0 = numpy.ones(nx)*background
    rho0[10:20] = 50
    return rho0



In [24]:

    
def integrate(rho0_background,t_f,rho0):
    """
    Integrate the system from time t0 = 0 to time t_f.
    
    Inputs:
    --- rho0_background, "background" in initial data function
    --- t_f, the final time we integrate to.
    --- rho0,  the boundary condition, rho(x=0) for all t
    """
    t = numpy.arange(0,t_f,dt)
    rho = numpy.empty((len(t),nx))
    rho[0] = initial_data(rho0_background)
    for i in range(1,len(t)):
        rho[i] = time_step(rho[i-1],rho0)
    return t,rho



In [25]:

    
def make_animation(x,t,rho):
    """
    Animate our traffic solution given an analytic solution
    rho(x,t) given by the arrays x,t,rho
    """

And integrate!



In [26]:

    
t,rho=integrate(BG_VAL,hours(7),BG_VAL)

Find indices for times we care about:



In [27]:

    
i3 = numpy.where(minutes(t)==3)[0][0]
i6 = numpy.where(minutes(t)==6)[0][0]
print i3,i6

Find velocities we care about.



In [30]:

    
print "Minimum velocity at time zero is: {} m/s".format(mps(numpy.min(V(rho[0]))))
print "Average velocity at time t=3 minutes is: {} m/s".format(mps(numpy.mean(V(rho[i3]))))
print "Minimum velocity at time t=3 minutes is: {} m/s".format(mps(numpy.min(V(rho[i3]))))
print "Minimum velocity at time t=6 minutes is: {} m/s".format(mps(numpy.min(V(rho[i6]))))









    



Minimum velocity at time zero is: 30.2222222222 m/s
Average velocity at time t=3 minutes is: 33.872218191 m/s
Minimum velocity at time t=3 minutes is: 30.9864026806 m/s
Minimum velocity at time t=6 minutes is: 34.4926790592 m/s

Make an animation for good measure



In [29]:

    
fig = plt.figure(figsize=(11,8))
ax = plt.axes(xlim=(0,L),ylim=(0,1.1*numpy.max(rho)))
line = ax.plot([],[],lw=2)[0]
time_text = ax.text(0.1*x[-1],0.9*numpy.max(rho),'',
                    transform=ax.transAxes)
def init():
    line.set_data([],[])
    time_text.set_text('')
    return line,time_text
    
def animate(i):
    line.set_data(x,rho[i])
    time_text.set_text('t = {}'.format(t[i]))
    return line,time_text
    
animation.FuncAnimation(fig,animate,frames=len(t),
                        init_func=init,interval=100)



In [29]: